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^ ; ABSTRACT 

. The rapid pace of extrasolar planet discovery and characterization is legitimizing the study 

of their atmospheres via three-dimensional numerical simulations. The complexity of atmo- 
spheric modelling and its inherent non-linearity, together with the limited amount of data 
Qh ' available, motivate model intercomparisons and benchmark tests. In the geophysical commu- 

P^] . nity, the Held-Suarez test is a standard benchmark for comparing dynamical core simulations 

pH ' of the Earth's atmosphere with different solvers, based on statistically-averaged flow quan- 

ta !■ tities. In the present study, we perform analogues of the Held-Suarez test for tidally-locked 

' ' exoplanets with the GFDL-Princeton Flexible Modeling System (FMS) by subjecting both 

^ . the spectral and finite difference dynamical cores to a suite of tests, including the standard 

' benchmark for Earth, a hypothetical tidally-locked Earth, a "shallow" hot Jupiter model and 

d . a "deep" model of HD 209458b. We find qualitative and quantitative agreement between the 

solvers for the Earth, tidally-locked Earth and shallow hot Jupiter benchmarks, but the agree- 
. ment is less than satisfactory for the deep model of HD 209458b. Further investigation reveals 

£> ' that closer agreement may be attained by arbitrarily adjusting the values of the horizontal 

dissipation parameters in the two solvers, but it remains the case that the iricignitucle of the 
horizontal dissipation is not easily specified from first principles. Irrespective of radiative 
transfer or chemical composition considerations, our study points to limitations in our ability 
to accurately model hot Jupiter atmospheres with meteorological solvers at the level of ten 
percent for the temperature field and several tens of percent for the velocity field. Direct wind 
C~ > ) measurements should thus be particularly constraining for the models. Our suite of benchmark 

tests also provides a reference point for researchers wishing to adapt their codes to study the 
atmospheric circulation regimes of tidally-locked Earths/Neptunes/Jupiters. 

^ \ Key words: planets and satellites: atmospheres - methods: numerical 

.5?; 

1 INTRODUCTION 

The nascent field of extrasolar planets is rapidly expanding, as evidenced by the flood of discoveries made in the past decade alone 
(e.g., Udrv & Santos) 2007 : Seager & Demingll2010[) . The abundance of exoplanetary data has legitimized several new fields of inquiry . 
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Amon g these is the theoretical study of exoplanetary atmospheres via detailed numerical simulations (for reviews see Showman et al„ 2008, 



2010), which describe the atmospheric dynamics, its radiative transfer, as well as — in principle — the chemistry and the cloud physics 

l, — „ jl„„j uz — — nu^Xj u„„ n i rin iu™ a u„„j r — — ttt — , iu^j rrr rr; — , jl„„j 



I Showman & Guillotll2002l: Icho et al lboqi 120081: Tcooper & Showmanll2005l bOQfj: iLangton & Laughlinll2008l: fMenou & Rauschetl !2009: 



Showman et al. 2009; Burrows et al. 2010; Rauscher & Menou 2010; Thrastarson & Cho 2010; Dobbs-Dixon et al. 2010). The complexity 
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Table 1. Table of Parameters and Commonly Used Symbols 



Quantity (Units) 


Description 


Earth (Held-Suarez) 


Earth (Menou-Rauscher) 


Shallow Hot Jupiter (Menou-Rauscher) 


HD 209458b 


Resolution 


conventional shorthand 


T63/G72 


T63/G72* 


T63/G72* 


T63/G72* 


N v 


vertical resolution 


20 


20* 


20* 


33* 


At (s) 


computational time step 


1 200 


1200* 


120* 


120* 


t' 1 (s _1 ) 


^ ■ j- - . + 


1.15741 x 10 — 4 


1.15741 x 10 — 4 


0.33422537 


0.32785918 




hyperviscous dissipation time 






10 ^ hot Jupiter day 


in — ^ vi n Tnoa^sti H'm 


K 


horizontal mixing coefficient* 


0.35 


0.35 


0.35 


0.1-1 


e 


longitude 


0-360° 


0-360° 


0-360° 


0-360° 


* 


latitude 


-90° -90° 


-90° -90° 


-90°-90° 


-90°-90° 


P s (bar) 


mean surface pressure 


1 


1 


1 


220 


Tfric ( da y) 


Rayleigh friction time 


1 


1 


oo 


oo 




planetary boundary layer 


0.7 


0.7* 






Trad ( da y) 


Newtonian relaxation time 


4-40 


15 


Tv/Slp sa 1.731 


equation IB 11 


Ti„it (K) 


initial temperature 


264* 


264* 


1800* 


1759 



Zstra (m) height of tropopause 1.2 x 10 4 2 x 10 6 $ 

a s t ra location of tropopause — ^;0.22 0.12 % 

Tsurf (K) surface temperature at equator 315 288 1600 % 

Tatra ( K ) stratospheric temperature 200 212 1210 % 

AT E p (K) equator-to-pole temperature difference 60 60 300 % 

AT B tra (K) tropopause temperature increment — 2 10 % 

AT Z (K) stability parameter 10 % 



Cp(Jkg _i K _i ) specific heat capacity at constant pressure 1004.64 1004.64 13226.5 14308.4 

7i(JItg _1 K _1 ) ideal gas constant 287.04 287.04 3779 4593 

k = TZ/c p — 2/7 2/7 2/7 0.321 

S2 p (s _1 ) planetary rotation rate 7.292 x 10 -5 7.292 x 10 -5 2.1 x 10~ 5 2.06 X 10" 

ffp(ms -2 ) planetary surface gravity 9.80 9.80 8 9.42 

Rp (m) planetary radius 6.371 X 10 6 6.371 X 10 6 10 8 9.44 X 10 7 



Note: unless otherwise stated, "day" refers to an Earth day (86400 seconds), 
j: Spectral models only. *: Finite difference models only. 
4b Value(s) used is different from in original publication. 
•V Value not explicitly specified in original publication. 
<0> : Vertical levels are logarithmically spaced. 
t : Thermal forcing of HD 209458b is given by equation (26]. 



of atmospheric modelling and its inherent non-linearity motivate clea n comparisons between studies that either utilize different methods of 

solution or even implement the same methods differently (Held 2005D . 

In the geophysical fluid dynamics community, it has long been suggested bv lHeld & Suarea ( 19941) that a useful comparison is between 
"dynamical cores", which are codes that deal with the essential dynamics of an atmosphere and omit details such as radiative transfer. In 
a pair of benchmark calculations. iHeld & Suarezl 1 1994) demonstrated that Earth-like simulations using two different methods — spectral 
and finite difference solvers — produced quantitatively similar profiles of the tempora lly-averaged, zonal -mean t emperature and zonal wind , 
as functions of vertical height (or pressure). Attempts have been made (e.g., Rauscher & Menou 20ld versus Cooper & Showman] 20061) 
to compare hot Jupiter models, but these were performed via different simulation platforms. In this study, we generalize Held-Suarez-type 
benchmarks to include tidally-locked exoplanets, using a single simulation platform — namely the dynamical cores of the Flexible Modeling 
System (FMS) . 

Our main finding is that while we can achieve qualitative and quantitative agreement for the Earth (and shallow hot Jupiter) tests, 
noticeable differences appear when we simulate the deep atmospheric circulation of the hot Jupiter HD 209458b. Closer agreement may 
be attained by specifying arbitrary values for the horizontal dissipation parameters — by trial and error — but it remains the case that 
the magnitude of the horizontal dissipation cannot be rigorously specified. Dynamical uncertainties at the level of > 10% therefore exist 
both between simulations utilizing different methods of solutions and also within the same method of solution, which may ultimately have 
implications for studies attempting to match observed versus simulated atmospheres of extrasolar planets. 

Operationally, we implement both the spectral and finite difference cores of the F MS and subject them to a battery of tests, including 



the Held-Suarez benchmark for Earth ( Q4. lb . a hypothetical tidally-locked Earth ( q|4.2 



Merlis & Schneider 



a " shallow" hot Jupiter 



mode l ( qT3llMenou & Rauschej|2009l) and a "deep" model for HD 209458b ( &3!4l ICooper & Showmanll2005L 12006, ; iRauscher & Menovj 
In Sj2] we discuss the governing equations handled by meteorological solvers such as the FMS. In ij3] we briefly describe the FMS. Our 
results are collectively stated in i|4]and we discuss their implications in fj5] Table[T]lists the parameters and commonly used symbols in our 
study, while Tablef2]describes the resolutions of the simulations. Appendices rAllBl andlClcontain technical details and useful fitting functions 
relevant to simulating the atmospheric circulation on the hot Jupiter HD 209458b. 



2 THE PRIMITIVE EQUATIONS OF METEOROLOGY 

The stud y of (terrestrial) meteo rology involves solving the Navier-Stokes and thermodynamic equations on a rotating sphere (e.g., 
Chapter 14 of Kundu & Cohenl 2004 ). Such an endeavour is usually inefficient or even intractable without invoking some simplifications, 



© 201 1 RAS, MNRAS 000. [711231 



Atmospheric circulation of tidally -locked exoplanets 3 



which results in a set of equa tions known as the "prim itive^] equations of meteorology" (e.g., ISmagorinskvlll963l . ll964l : Chapter 3 of 
Washington & ParkinsonlEooi : Chapter 2 of lVallisl f2006). The first simplification involves the assumption of vertical hydrostatic equilib- 



dP 



= ~99 



= —TIT, 



(1) 



dz ra dlnP 

where P denotes the pressure, z is the vertical/radial coordinate, p is the mass density of the fluid, g is the acceleration due to the gravity 
of the planet, <f> = gz is the geopotential, 1Z is the ideal gas constant and T is the temperature. The hydrostatic approximation filters out 
vertically propagating sound waves, but allows for vertically propagating gravity and Rossby waves as well as horizontally propagating waves 
in general. On large scales, hydrostatic equilibrium is a good approximation because the vertical pressure scale height H is much less than 
the planetary radius R p , 



rhgRp 



6 x 10" 



f_T\ /ft 



R* 



(2) 



V 1000 Kj \ 2m H 10 m s-2 R, 

where &b is the Boltzmann constant, m is the mean molecular mass, mn is the mass of a hydrogen atom and Rj « 71492 km is the 
(mean) radius of Jupiter. Such an ass umption precludes the explicit treatme nt of small-scale, three-dimensional turbulence, which may be a 
non-negligible source of dissipation jGoodmarJliooj iLi & Goodmalboioh . 

Consider the quantity r — R p + z. The second approximation then replaces r with R p in the equations of motion except where the 
former is used as the differentiating argument. The third approximation neglects the Coriolis terms in the horizontal momentum equation 
involving the vertical velocity. These approximations are collectively made such that angular momentum and energy conservation are ensured 
JVallisl2006t) . 

Let ve and denote the zonal (east-west) and meridional (north-south) components of the flow, respectively. The equations of mo- 
mentum and mass conservation are 



Dve 
Dt 

Dv<s> 
Dt 




2fiu$ sin $ + 



= — 2Qve sin $ 



«ev* tan $ 
Rp 
v% tan <J? 



dP 



Rn 



pR p cos $ 96 ' 
1 dP 

Jr~ p M' 



(3) 



(DP\ 



V.w = 0, 



dP V Dt J 

where v denotes the velocity vector. In a departure from traditional notation, we denote the latitude and longitude by $ and S, respectively. 
Equations Q3 and ([3} are augmented by the first law of thermodynamics, 

DT kT DP 



Dt 



P Dt 



+ — , 



(4) 



where k = lZ/c p and c p denotes the specific heat capacity at constant pr essure. The diaba tic heating is denoted by Q. For an ideal gas, 
c p — c v + 1Z, where c v is the specific heat capacity at constant volume. Goodmanl 1 2009b has remarked that equations {T}, <[3j and © 
collectively describe a frictionless heat engine, where no viscous terms exist to convert mechanical energy back into heat. 

Solving the equations explicitly with z is computationally awkward, especially when dealing with non-uniform topography. Instead, P 
is used in place of z such that the temporal derivative following the flow is 

vq d v<s> d DP d 
! ~DtdP' 



D_ _ d_ 

Dt ~ dt + R v cos $ dO 



+ Rp"d$ + 



(5) 



(6) 



In addition, the pressure is normalized by the instantaneous surface pressure P 3 , such that 

P 

This is also known as Phillips' cr-coordinate and was designed to deal with mountainous terrain in geophysical calculations JPhillipsll 19571) 
By definition, the a = 1 level tracks the (exo)planet's orography (if any). 



3 THE GFDL-PRINCETON FLEXIBLE MODELING SYSTEM 

The Flexible Modeling System (FMS) is an open source, parallel simulation platform developed at the Geophysical Fluid Dynamics 
Laboratory (GFDL) of Princeton University. The FMS has three core options: finite difference, spectral and finite volume. In this study, 
we implement the Memphis release of FMS and utilize both the spectral and finite difference ("B-grid"|^| dynamical cores. As the FMS 
utilizes MKS units, some of the discussion in the paper will follow suit. In this section, we describe some salient features of the FMS. 



1 Fro m a historical viewpoint, the term "primitive" is a misnomer, since it means "full" rather than "simple" (see Chapter 3.2 of Washington & P arkinson! 
l2005h , 

2 We note t hat the B-grid finite di f ference scheme is an old one and is k nown to be less accura te than the more commonly-used "C-grid" scheme, e.g., as 
employed bv lHeld & Suarej h994l) . ICooper & Showmanl J200ll2006l) and lShowman et all <2009h . 
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Readers interested in more technical details may consult http : / / www . qf dl ■ noaa . qov/f ms.F or an overview of the various simulation 
platforms which are publicly available, please refer to Chapter 5 of Washington & Parkinson! 1 20051) . 



3.1 Spectral Core 

In th e spectralfl dynamical c ore of the FMS, the hydrodynamic variables are described as a sum of spherical harmonics truncated at 
Nh terms ( Gor don & Sternlll982h . Triangular truncation is used in our implementation of the FMS, such that the truncation is rotationally 
symmetric — a function and its rotated counterpart are both expressible within this truncation (see §13.6.2 of Holtor] 2004 ). On a sphere, 
the number of zonal and meridional waves retained are Nh and Nh + 1, respectively, to prevent aliasing. The corresponding number of 
longitudinal grid points (7Vi on ) is always twice that of the latitudinal grid points (JVi a t). Domain decomposition is ID in the spectral core: the 
number of processors allocatable to computing a given model is A r i at /2. 

A key aspect of any spectral model is spectral blocking, which is the accumulation of n umerical noise — specifically, enstrophy — at the 
smallest grid scales, since spectral codes are intrinsically non-dissipative ( Stephenson 1994). Numeric al "hyperviscos ity" is needed to mimick 
enstrophy dissipation at the smallest length scales, analogous to a two-dimensional turbulent cascade Jshapirclll97lh . The hyperviscous term 
takes the form, 



2> 



hyper 



(V x v) z 



(7) 



where n v i s the hyperviscosity dampi ng order and (V x v) z is the relative vorticity. Following [Held & Suarea dl994T ). lMenou & Rauscher 
(2009) and lRauscher & Men ou (2010), we adopt n„ = 4. Within the FMS, one may either specify the hyperviscosity coefficient (v) or the 
dissipation rate (~ vV 2n ")\ we will discuss this issue further in ij3.3l The spectral core has an optional switch to ensure global energy 
conservation, which we enforce for all of our simulations. 

The conventional shorthand notation used to describe the resolution of the spectral models is T7VhLiV v , where Nh is the horizontal 
resolution while jV v is the number of vertical levels. The fiducial resolution we will adopt for our spectral simulations is T63, which cor- 
responds to AW = 192 an d N\^t = 96. Finite differencing is used for the vertical grid by employing the Simmons-Burridge scheme 
(Sim mons & Burridgelll981f) . For example, the lowest layer modelled within such a scheme has a = 0.95-1; the boundaries between the 
layer are called the "half levels" (i.e., a = 0.95 and 1). A noteworthy feature of the Simmons-Burridge scheme is that the simulation output 
is not exactly at the midpoint between the half levels (i.e., not at a = 0.975 in this example). Therefore, it should be noted that when we 
present our results, we usually quote the pressure level P as the larger of the pair of half level values (e.g., for P — 0.95-1 bar layer, we 
label it "P = 1 bar"). 

A key advantage of the spectral method described here is that it does not require special (damping) treatment at the poles. This is not 
the case for the finite difference core. 



3.2 Finite Difference Core 



The finite difference dynamical core of the FMS uses a "Arakawa B-grid" (see Chapter 4.2 of Washington & Parkinson 2005) for the 
horizontal coordinates, which belong s to a family of finite difference grids where the temperature and velocity are solved at staggered points 
l Wvmanll996l : I Anderson et alj2004 ). The vertical grid uses a hybrid a-P coordinate system; the labelling of the different model layers again 



follows the larger of the pair of half level values (see i)3.U . Finite differenc ing is used for both the horizontal and vertical grids. Analogous 
to the case of the spectral core, small-sca le noise accumulates in the B-grid ( Sha piro! 197fj) and has to be damped via a "horizontal mixing'Q 
algorithm dRoeckner & von Storchlll980|) . In the finite difference core, the second-order operator for horizontal mixing is defined as 



« = -x^Kp. E ^ [ Ki % (J 7 )] , (8) 

where T denotes the temperature or zonal/meridional velocity components (i.e., T or v), Ajr is the area of each grid box (for either the 
temperature or velocity fields), 

Ki = K A, Ji{Ar) Ji(APj), (9) 

and the index i = $, O. Denoting an arbitrary quantity by Q, the operators Z»(Q) and Ji(Q) yield the difference and average between 
adjacent grid points along the i-axis, respectively. The difference in pressure between half-levels at an index j is 

AP 3 ■ =P J+1/2 ~P^ 1/2 . (10) 

The quantity A.; is a constant that describes the strength of the horizontal mixing with latitude and must satisfy the numerical stability 
condition: A 4 /C ^ 1/8. It is important to note that the operator defined in equation {8]l is essentially a Laplacian and is applied twice to the 



3 Strictly speaking, the code is pseudo-spectral because only the linear terms in the governing equations are transformed to the spectral domain, while the 
non-linear terms are computed on a finite difference grid. This statement is independent of the method of solution for the vertical coordinate. 

4 Also termed "horizontal diffusion". 
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-H n 



(11) 



temperature and velocity fields, 

dT _ 
dt ~ At 

where At denotes the time step, implying that the horizontal mixing scheme is fourth order in nature^ The horizontal mixing coefficient has 
a range of values of K ^ 1. Its default value within the FMS is AC = 0.35, which we will adopt throughout unless otherwise stated. 
Damping is increased for |$| > 80° but is uniform with longitude. 

The common problem faced by any finite difference code which solves the fluid equations on a sphere using the longitude-latitude 
coordinate system is that, for numerical convergence to be attained, the minimum time step needed is proportional to the zonal grid spac- 
ing A 8 via the Courant-Fredricks-Levy (CFL) condition. This implies that At — > towards the poles. For example, iDobbs-Dixon et al 
(2010) truncate their latitudinal grid at $ = ±70° in their simulations. To alleviate this problem, a technique known as "polar filtering" is 
applied at high latitudes to damp the shortest resolvable waves such that a non-zero time step can be taken I Shapiro 1971 ; Asselinl 1972 ; 
iTakacs & Balgovindl ll983). Like in the spectral core, the finite difference core has a switch to ensure global energy conservation, which we 
set to "on" for all of our simulations^ 

The shorthand notation used for resolution is GTVhLiVv where Nh now refers to half of the number of grid points in longitude (i.e., 
around a latitude circle). Alternatively, one can use the notation NN-^LNv where is the number of latitudinal points between the 
north/south pole and equator. Domain decomposition is 2D in the finite difference core. The fiducial resolution for our finite difference 
simulations is G72/N45, which corresponds to iVi on = 144 and A r i at =90. 

Finally, we note that an alternative approach within grid-based meth ods is to adopt a "cubed- sphere" grid, w hich circumvents the 
problems at the poles at the price of dealing with a non-orthogonal grid (e.g., Adcroft et alJl2004 : Showman et al.112009}) . 



3.3 Horizontal Dissipation 



It is important to note that the horizontal dissipation schemes described above are reasonably well-motivated but nevertheless non- 
rigourous. Hyperviscosity and horizontal mixing are numerical tools unsupported by any fundamental physical theory, yet are routinely used 
by research groups studying terrestrial and exoplanetary atmospheric circulation. There is no rigourous way to choose their magnitudes. On 
Earth, the magnitude of horizontal dissipation can be calibrated on the basis of the known flow, but this is not (yet) — and may never be — 
the case for extrasolar planets. The use of horizontal dissipation is related to the notion that turbulent cascades of hydrodynamical quantities 
(energy, enstrop hy) are only parti ally modelled. As such, both v and K, should be regarded as free parameters in any model of atmospheric 
circulation (e.g., ISteph enson 1994). 

In practice, the pragmatic aim is to dissipate small-scale numerical noise within a fraction of a planetary rotation (i.e., one day on a 
tidally-locked planet). For the spectral core, the dissipation time on the scale of a resolution element is 

where n„ = 4 is usually adopted. Spectral simulations with smaller values of t v are generally more dissipative. To meaningfully compare 
spectral simulations with different numerical resolutions, we need to keep the hyperviscosity v fixed by using equation d!2t to scale the 
dissipation rat^| assumed, 



For example, the dissipation rate used in the T63L33 run, for the deep model of HD 209458b, is t~^63 ~ °-33 s" 1 (» 10" 5 HD 209458b 
day). Therefore, a T31L33 run would use t v ~ 3 x 10 -3 HD 209458b day. Within the FMS, specifying the dissipation rate as an input 
parameter is thus termed a "resolution dependent" run. Alternatively, specifying the hyperviscosity v constitutes a "resolution independent" 
run. In general, we find that dealing with a dissipation rate (with units of s -1 ) is somewhat more intuitive than having to vary v (with units 

of m 8 s" 1 ). 



In the original lHeld & Suarezj I ] 19941) T63L20 (JViat = 96) spectral simulations, the dissipation rate was chosen to be 1.15741 x 10 4 
s -1 w hich corresponds to a dissipation time of about 0. 1 Earth day. For their T42L15 (iVjat = 64) Earth-like simulations. lMenou & Rauscher 
l iooi) use v = 1 . 18 x 10 37 m 8 s 1 , which corresponds to t v ~ 9 x 10 3 Earth day. For their T42L15 hot Jupiter simulations, 
Menou & Rauschei 1 20091) u se v = 6.28 x 10 47 m 8 s -1 , which corresponds to t v 
simulations of HD 209458b. iRauscher & Menoul J20ld) use v = 8.54 x 10 47 m 



day. For our Earth-like simulations, we choose t v = 0.1 day following iHeld & Suarezl ( 11994b . For our hot Jupiter simulations, we choose 



- 2 x 10~ 4 hot Jupiter day. For their T31L33 (AT lat = 48) 
s" 1 which i s equivalent to t v ~ 9 x 10 -4 HD 209458b 



5 The Memphis release of the FMS uses defaults of second and fourth order for the wind and temperature horizontal mixing schemes, respectively. We have 
performed two separate suites of simulations where the wind scheme is set to second or fourth order and find little difference between the Held-Suarez statistics 
generated. 

6 We note that the application of horizontal mixing and polar filtering result in small violations to the conservation of mass and energy. 

7 From an operational standpoint, we note that when the dissipation rate is too small, the simulations will crash even when very small time steps (e.g., At = 1 
s) are taken. Therefore, there is a practical upper limit to the value of t„ assumed. 
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Table 2. Table of Simulation Resolutions 



Simulation 


Spatial Resolution 


Angular Resolution 


Examples (3D) 


T21 


64x32 


5.625° 


HMP 


T31 


96x48 


3.75° 


Burrows et al. (20101: Rauscher & Menou (2010); HMP 


T42 


128x64 


2.8125° 


Menou & Rauscher (2009); Thrastarson & Cho (2010)+ 


T63 


192x96 


1.875° 


Held & Suarez ( 1994); HMP 


G24 


48x30 


(7.5°, 6.0°) 


HMP 


G36 


72x45 


(5.0°, 4.0°) 


Cooper & Showman (2005. 2006); HMP 


G48 


96x60 


(3.75°, 3.0°) 




G72 


144x90 


(2.5°, 2.0°) 


Held & Suarez (1994); Dobbs-Dixon et al. (2010)*: HMP 



Note: The acro nym "HMP" refers to the p resent study. 

t: lMenou & Rauscheil J2009h and lThrastarson & Choi J201ph presented mainly T42 models 
but al so examined T85-T170 and T21-T85 ones, respectively, for convergence tests. 
|: Dobb s-Dixon et alj «2010f> used a resolution similar to G72 for their simulations 
(160 X 64; 2.25° X 2.1875°) but truncate their latitudinal grid at $ = ±70°. 



t v = 10~ 5 hot Jupiter day such that t o within a factor of a few, our chosen value for t v is consistent with those used bv lMenou & Rauscher 
1 20091) and Rauscher & Menou! J201dh . Table[JJlists our choices for t ~ 1 , which were made to match the values in the original publications as 
closely as possible, while bearing in mind that all of these choices have no strict justification beyond the requirement that the model can be 
integrated without the detrimental accumulation of small-scale noise. 

For the finite difference (B-grid) core, the horizontal mixing coefficient /C plays the analogous role of t v — and not v — in the spectral 
core (see Appendix[C]l- Its default value within the FMS is /C = 0.35, which we adopt unless otherwise stated. Varying K, is in effect changing 
the value of the analogue of t v . Conversely, keeping /C fixed and varying the resolution of the simulation effectively varies the value of 
the analogue of v. We are unable to write down a simple analytical expression relating K, and v, but note that it is possible to measure the 
analogue of t v in a finite difference simulation (see Appendix|Cj. While a correspondence between the dissipation parameters in the spectral 
and finite difference cores may exist, we do not consider it to be straightforward. Our main intention is to demonstrate that it is possible to find 
equivalent pairs of values for t v (or v) and IC by trial and error (see £|4.4b . which has implications for researchers wishing to adapt existing 
simulation platforms implementing different solution methods (and numerical dissipation schemes) to study the atmospheric circulation on 
exoplanets. 



3.4 Initial Conditions 

The default initialization in both the spectral and finite difference cores uses the simplest assumption: isothermality with no wind. Every 
temperature point on the solution grid is set to T = Ti n i t where Tj n jt is an initial temperature and may be regarded as a free parameter. 
The tests we will describe in S|4]use values of Tinit which are tabulated in Table [JJ We will see that the active (t~j 7^ 0) layers of the 
atmosphere, where the temperature rapidly relaxes towards T oq , are somewhat insensitive to the choice of Tinit and produce results that are 
broa dly consistent with previou s studies. The simulations are started with a small initial perturbation in the vorticity field. 



lly consistent with previous studies, the simulations are started with a small initial perturbation in the vorticity held. 
Thrastarson & ChoU201ol) have argued that initializing the simulations with non-zero winds can produce both qualitative and quantita- 



tive differences in the results, because the applied thermal forcing is projected differently onto the normal modes of the atmosphere under 
different initial wind conditions. We consider this issue to be beyond the scope of the present study. 



4 ATMOSPHERIC DYNAMICAL CORES: TESTS FOR EARTH AND HOT JUPITERS 
4.1 Held-Suarez Benchmark Test 



As a first check on our computational setup, we reproduced the iHeld & Suarea dl994h benchmark test with both the spectral and 
finite difference versions of the dynamical core. The effects of stellar irradiation, geometry, etc — known as the "thermal forcing" — are 
encapsulated in the "equilibrium temperature" function, 



T cq = max{T stra ,T H s}. 



where T stra = 200 K is the stratospheric temperature, 



T HS 



T surf - AT E p sin' $ 



AT Z In ( — ) cos' <E> 



(14) 



(15) 
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Figure 1. Zonal-mean temperature, temporally averaged over 1000 days, for the Held-Suarez benchmark test for Earth. Contour levels are in units of K; 
P s = 1 bar. Left: T63L20 spectral model. Right: G72L20 finite difference model. 

T63L20: zonal — mean zonal wind G72L20: zonal-mean zonal wind 




$ (deg) 4> (deg) 

Figure 2. Zonal-mean zonal wind, temporally averaged over 1000 days, for the Held-Suarez benchmark test for Earth. Contour levels are in units of m s _1 ; 
P s = 1 bar. Left: T63L20 spectral model. Right: G72L20 finite difference model. 



Ts Ur f = 315 K is the surface temperature at the equator and ATep = 60 K is the equator-to-pole temperature difference. The parameters in 
the preceding equation specific to the Held-Suarez forcing are set to be 

AT Z = 10 K, 

(16) 

Po = 1 bar, 



while the other parameters of the test are described in Table Q] The initial temperature is not specified in iHeld & Suarea i 1994) . but the 
default value in the spectral code is Tmit = 264 K; we will adopt this value for both the spectral and finite difference simulations. 
The FMS implements a simple Newtonian relaxation of the temperature field, where a damping coefficient, 

1 I 0, a < a b , 

^Newton = h < / , , \ / „ „ \ . (17) 

is applied to the temperature field relative to equilibrium (T cq — T). In the original implementation of Held-Suarez forcing, we have r ra d,u = 4 
days, r ra d,d = 40 days and at = 0.7 denoting the top of the planetary boundary layer in <r-coordinates. Later, we will also implement only 
a single value of the Newtonian relaxation time, i.e., r ra d = r ra d,u = T" ra d,d, such that "Newtonian cooling" is represented by the term 

T cq — T 

QNewton = • (18) 

^~r a d 

Low-level winds are damped on a time scale rfric using the damping coefficient, 

I 0, a ^ at, 

V Ra . ylelgh = \ ' d9) 

Such a prescription is known as "Rayleigh friction" (or drag) and mimicks boundary-layer friction between the atmosphere and the surface 
of the Earth. In the FMS, Rayleigh friction is applied to the velocity field: — ©Rayioigh^- Note that <7{, 7^ <7 str a in general, where cr stra is the 
location of the tropopause, the transition layer between the troposphere and stratosphere. We will later implement thermal forcings which are 
different from equation j 14b . 
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For the spectral model, the resolution used for the published results of Held & Suarezl 1 1994 ) is T63L20. The corresponding longitude 
versus latitude grid for this resolution is 192 x 96, which allows the model to be computed on up to 48 processors simultaneously. The 
default setting for the FMS Held-Suarez benchmark uses a hyperviscous dissipation rate of 1.15741 x 10~ 4 s _1 (~ 0.1 day). The normalized 
pressure (0 < a ^ 1) is equally spaced with 20 vertical levels. 

For the finite-difference model, the resolution used in the published results o f lHeld & Suarezl 



1994) is G72L20 (144 x 90). The vertical 



levels are treated with a hybrid a-P coordinate system: the terrain-following a-coordinate is used near the planetary surface and transitions 
to the P-coordinate well above the surface. The default setting for the horizontal mixing coefficient is K, — 0.35. 

Figures[TJand[2]show the zonally-averaged (or zonal-mean) temperature (in K) and zonal wind speed (in m s _1 ), respectively. Following 
Held & Suarezl d 19941) . we ran both sets of simulations for ttotal = 1200 days, but discarded the first tdiscard = 200 days in order to 



eliminate features due to the different initialization schemes. We see that our results are consistent with those presented in Figures 1 and 2 of 



Held & Sua rez ( 1994); the spectral and finite difference results are also in broad quantitative agreement. For the rest of the paper, we refer to 



the quantities shown in Figures\T]and\2\( i.e. , 1000-Earth-day averages of zonal-mean profiles) as the "Held-Suarez statistics". 
We conclude that our implementation of the Held-Suarez benchmark test for Earth is successful. 



4.2 Hypothetical Tidally-Locked Earth Benchmark Test 

Thermal forcing for a tidally-locked (exo )planet can be mimicked by replacing the — sin 2 $ term in equation (|15[l with a term that is pro- 
portional to + cos(9 - 180°) cos $ (see, e.g. JCooper & Showmanl2005l ; lMenou & Rauscherl2009l:lMerlis & Schneiderl20ld) . Additionally, 
for a hypothetical tidally-locked Earth the rotation rate has to be reduced, 



- ^ n 



fip/365, 



(20) 



such that one planetary day is equal to one planetary year. As a prelude to simulating the atmospheric circul ation on (tidally-locked) ho t 
Jupiters, we first examine the case of such a tidally-locked Earth at 1 AU. Such a case study was conducted bv lMerlis & Schneider 
who considered more sophisticated physics than is the case for our dynamical core simulations, including an active hydrological cycle, a gray 
radiative transfer scheme with a pressure-dependent opacity, and an explicit formulation for atmosphere-surface exchanges on an aquaplanet. 
We implement the following thermal forcing 



T oq = max {T st 
Ths = 



,Ths} ■ 



T surf + AT EP cos (6 - 180°) cos $ - AT Z In ( — 



(21) 



placing the substellar point is at (O — 180°, $ = 0). We adopt the same set of parameters as described in Table [JJ for the Held-Suarez 
benchmark, including for the implementation of Newtonian relaxation and Rayleigh friction. The time step used is At — 600 s. 

The first row of Figure [3] s hows snapshots of t he temperature field at Day 1200 and a = 1.0. The temperature field, which should be 
compared to Figure 1 o flMerlis & Schneidej hold) , shows an atmospheric temperature structure dominated by radiative forcing rather than 
advection. The second, third and fourth rows of Figure [3] show the temporally averaged (over 1000 days) zonal wind profiles, as functions of 
longitude and latitude, at o = .25, 0.55 and 1.0, respectiv ely. These values of a were chosen to match as closely as possible the a — 0.28, 
0.54 and 1 .0 values adopted by Merlis & Sdmeide^j20ic|) in the lef t column of their Figure 4. Despite our much simpler setup, our results 
are in qualitative agreement with those of lMerlis & Schneider! 1 2O10l) . showing the presence of a large, direct circulation cell centered on the 
substellar point. Furthermore, our spectral and finite difference simulations are in general agreement with a clear indication that discrepancies 
start cropping up near the poles, as may be expected because of the difficulties in treating the poles in the finite difference core. 

Figure[4]s hows the temporally averag ed meridional wind profiles at a — 0.25, 0.55 and 1.0, and s hould be compared to the ri ght column 
of Figure 4 of iMerlis & Schneider 1 2010h . We again attain qualitative agreement with the results of Merlis & Schneider (2010). capturing 
the large circulation cell centered on the substellar point, which exhibits poleward and equatorward motions in the dayside and nightside 
hemispheres, respectively, i.e., the atmosphere flows from the day to the night side. 

Our results in this subsection provide a useful prelude to the study of hot Jupiter atmospheres, because the simulation of a tidally-locked 
Earth is less computationally demanding. The benchmark tests described in this sub-section are thus an efficient way of check i ng on e's code 
before moving on to the hot Jupiter benchmarks. For variations on a theme of Earth-like models, please refer to lHeng & Vogt (2010). 



4.3 Benchmark Tests for Earth and "Shallow" Hot Jupiter 

Menou & Rauscher 1 20091) simulated the atmospheric circulation on Earth and hot Jupiters using the 3D IGCM spectral code with equal 



spacing in a (N v — 15). In the hot Jupiter case, this model is considered to be "shallow" because of the limited depth of the 3D atmosphere 
modelled (d own to only 1 bar). S uch a 3D model should not be confused with the " shallow water" l Menou et al. 2003b or "equiv alent 



barotropic" I Cho et al]|2003 , 200 8h models, both of which are essentially 2D. (See also Longuet-Higgins 
Heng & Spitkovskvll2009l .) 



1968, 



Kundu & Cohen 2004 and 
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Figure 3. Simulation of tidally-locked Earth. Top row: snapshot of the temperature (represented by colours) field at 1200 Earth days and cr = 1.0. The second, 
third and fourth rows are the temporally averaged zonal wind profiles at a = 0.25, 0.55 and 1.0, respectively. The left and right columns are for the spectral 
(T63L20) and finite difference (G72L20) simulations, respectively. Temperatures are in K and wind speeds are in m s . 
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Figure 4. Same as the second, third and fourth rows of Figure|3] but for the temporally averaged meridional wind profiles. 



4.3.1 Earth-like 

For an Earth-like setup, instead of equation J14t . lMenou & Rausc her (2009) employed 



1 2 

Tcq = Tvert + (Strop AT B P I - - Sm $ 



for the thermal forcing, where 



and 



l-op (Zstra + " g Stra ) + | ^ trop(^ btra^j _|_ \ 



1/2 



urf Ttrop^stra ~~h AX'stra, 



■2 ^ ^stra, 
•2 )!> £stra, 



->trop 



sm 
0, 



T^rfCT-CTatra) 1 
[ 2(l-CT s tra) J 



Z ^ Zstra Or a (Tstra, 
« > Z stra Or (7 < <X s tra- 



(22) 



(23) 



(24) 
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A single value of th e Newtonian relaxat i on tim e is considered (r ra( j 
in cr-coordinates. In 



15 days). We reiterate that <7 s tra is the location of the tropopause 
Menou & Rauschei t he in itial temperature is not explicitly specified, so we choose T^t = 264 K following 



Held & Suarel Jl994h . Also. llVlenou & Rauschej J2009t) apply Rayleigh friction only to the bottom-most layer of their T42L15 simulation 



where as we choose crj, = 0.7 for our Rayleigh friction scheme. An important difference between the two schemes is that lMenou & Rauschei 
i 2009h apply Rayleigh fr iction to the vorticity a nd divergence fields, while Rayleigh friction is i mplemented in the FMS as applying directly 



to the velocity field. iMenou & Rauschen J20091) consider this test to be a simplified version of the lHeld & Sua rez ( 1994) benchmark. 



Figure [5] shows the Held-Suarez statistics for this benchmark test. The label led contours are different from those in Figures Q] and 
[2] so as to facilitate direct comparison with Figure 2 of Menou & R auschei 20091) . Our results are temporally averaged over 1000 days, 
from day 200-1200, while IMenou & Rauschei J2009I) present results for day 150 The temperature profiles computed in our spectral and 
finite difference simulation s are essentially identical; they also match the temperature profile presented in the bottom panel of Figure 2 of 
Menou & Rauschei J2009I) . There are some noticeable differences between the zonal-mean zonal winds computed by our spectral and finite 



difference simulations, yet to a good degree they are mutually consistent and also agree with the top panel of Figure 2 of lMenou & Rauschei 
I 2009T) . 



4.3.2 Hot Jupiter 



i 2009h 



The shallow hot Jupiter model considers evenly spaced a levels where P s = 1 bar. The thermal forcing implemented bv lMenou & Rauschei 



T cq = T vcrt + /3 trop AT E p cos (6 - 180°) cos $, (25) 

placing the substellar poin t at (0 = 180°, $ = 0). No Rayleigh friction/drag is implemented, following Menou & Rauschei 1 20091) . Since 
Menou & Rauschei i2009) do not specify their choice of the initial temperature, we (arbitrarily) adopt Ti n i t = 1800 K. The Newtonian 
relaxation time is half a hot Jupiter day, ir/il p « 1.731 Earth days, which is short enough that the results appear insensitive to the choice of 
Tinit. Results produced with T- ln i t — 264 K and 1470 K are very similar to those reported below. 

Both the spectral and finite difference simulations of hot Jupiters were performed using a time step of At — 120 s, a factor of 10 smaller 
than for the (non-tidally-locked) Earth-like simulations. As before, we use ttotai = 1200 and fdiscard = 200 Earth days, both of which are 
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Figure 6. Results of spectral (left column; T63L20) and finite difference (right column; G72L20) simulations for the shallow model of hot Jupiters. Top row: 
temperature fields at 100 hot Jupiter days (ss 346 Earth days) and a = 0.7. Middle row: zonal-mean temperature, temporally averaged over 1000 Earth days. 
Bottom row: zonal-mean zonal wind, temporally averaged over 1000 Earth days. Temperatures are in K and wind speeds are in m s . 



sufficiently long for the simulation to reach quasi-equilibrium; we will retain these values for ttotal and fdiscard for simulations involving hot 
Jupiters throughout the paper. 

Figure [6] shows the usual Held-Suarez statistics, as well as snapshots of the temperature field at 100 hot Jupiter day s (« 346 Earth 
days) and a — 0.7. These snapshots are meant to be compared to the top panel of Figure 3 of iMenou & Rauschen d2009f) — we can see 
that there is general agreement with the qualitative features of the flow and the range of temperatures produced. It is unsurprising that the 
snapshots are not identical as the time required to reach quasi-equilibrium is slightly different for each simulation — therefore, one cannot 
attribute differences in the snapshots (top row of Figure [6} to the different methods of solution (spectral versus finite di fference). A more 
meaningful comparison is between the Held-Suarez statistics produced, both with Figure 2 of Menou & Rauschei] d2009h and between our 
pair of simulations (middle and bottom rows of Figure|SJ, where we witness general agreement. In general, this point should be kept in mind 
when examining simulation snapshots produced by different methods of solutions. 

Slight quantitative differences between the Held-Suarez statistics produced by the spectral and finite difference simulations exist. For 
example, the temporally averaged, zonal-mean zonal wind speed ranges from -630 m s _1 to 1239 m s~ x for the spectral simulation, but is 
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Figure 7. Newtonian relaxation time and day/night side temperature profiles as computed for HD 209458b. Left: diamonds represent the calculations from 
Figure 4 of Iro et al. (2005), while the dotted curve is our 4th order polynomial fit. Right: day (Td ay ) and night (T n i s ^ t ) side temperature profiles computed 
from the globally-averaged profile of Iro et al. (2005; 7i ro ). 



-554 m s" 1 to 1096 m s 1 in the finite difference simulation. We will not pursue the cause of these differences for the shallow hot Jupiter 
model, which are probably due to a combination of resolution and choice of the magnitude of the horizontal dissipation. Instead, we feel that 
the "deep" model of HD 209458b provides a more meaningful exploration of these issues, which we will examine in £14.41 



4.4 Deep Benchmark Test for Hot Jupiter: HD 209458b 

Rauscher & Menou 1 2O10h modelled atmospheric circulation on the hot Jupiter HD 209458b from P 



1 mbar down to P 



220 



bar with N v — 33 unevenly spaced vertical levels. They find an upper atmosphere dominated by radiative forcing, due to the short radia- 
tive time scales, and an a dvection-dominated lower atmosphere with a low level of varia bility. General features of the flow (Figure 1 of 
Rauscher & Menoull2010l) agree qualitativ e ly with Figu re 1 of Cooper & Showmai] 1 2005), who employed the finite difference (instead of 
the spectral) method. [Cooper & Showman] 1 2005 , 20061) also modelled the atmosphere of HD 209458b down to deeper levels: P = 3 kbar 

instead of 220 bar. 

Key differences in the results of ICooper & Showmanl 1 2005 . 20061) and Rauscher & Menou 1 20ld) are of interest because they may 
indicate limitations in the different methods of solutions for the hot Jupiter regime. We summarize the differences: 



(i) The super-rotating e quatorial jet descends down to onl y about 7 bar in the simulations of Rauscher & Menoi] 1 20ld) . but can be found 



down to 50 bar in those of 



r Cooper & Showmanl J2005l . l200d) : 
(ii) The simulations of Cooper & Showman (2006) show predominantly super-rotating and counter-rotating flows in the upper and lower 



atmospher e, respect i vely. B y contrast. lRauscher & Menou (2010) find flows in both directions throughout the active layers of the atmosphere. 
Rauscher & Menou! ( 120100 interpret this difference as being due to the deeper reservoir of inert (r ra( j = oo) atmospheric layers in the models 



of 



Cooper & Showm an (2006), which counter-balances the angular momentum of the super-rotating wind higher up in the atmosphere. 



Using both the spectral and finite difference cores of FMS, we will see that the discrepancies described above vanish, implying that 
they probably arise from differences in initial/boundary conditions as well as setup. However, our initial simulations reveal other quantitative 
differences, which we will discuss. 



4.4.1 Setup 



To perform simulations spanning several orders of magnitude in a or P, one needs to implement uneven vertical spacing, the technical 
details of which are described in Appendix [A] Figure lATl il lustrates the setup ne eded for a simulation with N v — 33 to cover 1 mbar < P ^ 
220 bar, similar to the simulation of HD 209458b bv lRauscher & Menoul ( hold) . 

J2005b 



The other ingredients needed are the functional forms of r ra d and T cq . In the case of HD 209458b 



Iro ct al 



the Newtonian relaxation time in their Figure 4 using ID, time-dependent, radiative transfer models. Rauscher & Menoul J2OI0I) apply this 
calculation of r ra( j when P < 10 bar. Appendix |B]contains a polynomial fit to r ra d = Trad (J 3 )- The specification of the radiative relaxation 



time in turn specifies the "active" (T rad > 0) and "inert" (r 



0) layers of the atmosphere. 



The temperature profile for the thermal forcing of HD 209458b is given by equation (2) of Cooper & Showmanl ( 2005), 



[r n i g ht 



Tn'ight) C0S(9 



180°) cos$ 



1/4 



T 
- 1 1 



night 1 



90° < $ < 270° 
otherwise, 



(26) 



where T n i g ht and Tday are the temperature profiles as functions of pressure on the night and day sides, respectively, llro et al.l (12005) have 



© 2011 RAS, MNRAS 000,[Hl23] 



14 Heng, Menou & Phillipps 



396 499 603 706 810 913 1017 1120 1224 1327 1431 1075 1139 1203 1267 1332 1396 1460 1525 1589 1653 1718 




(deg) e (deg) 



1545 1597 1650 1702 1755 1808 1860 1913 1965 2018 2071 1698 1724 1751 1777 1804 1830 1857 1883 1910 1936 1963 




8 (deg) (deg) 



Figure 8. Snapshots of the flow field at about 340 HD 209458b days at P = 2.13 mbar (top left), 216 mbar (top right), 4.69 bar (bottom left) and 21.9 bar 
(bottom right). The T31L33 simulation was performed using the FMS spectral core with i„ = 10 — 5 HD 209458b day. Colors indicate temperature in K, while 
arrows represent the velocity vectors. 



computed the globally-averaged (between night and day ) temperature profile ( Tiro; s ee Appendix |Bj, which Rau scher & Menou 1 20101) have 
used to calcula te Tnight and Td ay ; see also Figure 1 of ICooper & Showman! ( I200q). We recompute these pro files using the calculations of 
Iro et al.l 1 120050 by solving the transcendental equation for T n i g ht (equation [22] of lCooper & Showmanll200q) . 

42£o = 3T 4 gM + (T night + AT cq ) 4 , (27) 

at each value of P. The temperature difference between the night and day sides, AT oq , is set equal to 1000 K for P ^ 1 mbar and 530 K at 
P — 10 bar. In between, AT oq is equally spaced at uniform intervals in log P. Our polynomial fits for T n i g ht and Td ay are given in Appendix 

m 

FigurefTJshows the radiative relaxation time and thermal forcing function used in our simulations of HD 209458b. The initial temperature 
is se t to Ti TO (P = lObar) = 1 759 K, which is the value of T cq at P = 10 bar. Our use of a constant initial temperature is simpler than 
what Rauscher & Menou 1 2010h implement, which is Tinit = Tnight for P < 10 bar and Ti n i t = Ti ro otherwise. In their T31L33 simulation, 



Menou (z 
erl J2009n~ 



Menou & Rauscher! d2009l) used 7Vi at = 48 such that small-scale numerical noise is dissipated on a time scale t v ~ 9 x 10 -4 HD 209458b 
day. For the spectr al simulations, we se t the d issipation rate to be exactly t^ 1 — 0.32785918 s _1 (t v = 10 -5 HD 209458b day) which is 
Rauscher & M enou feoioh value (to within a factor of a few) if the scaling relation in equation J 1 3 1 > is considered. For 



consistent with the 

the finite difference simulations, we initially adopt /C = 0.35 for the horizontal mixing coefficient. We will explore variations in t u and K, as 
well as in and N v . 



4.4.2 Results 



Figure [8] shows snapshots of the temperature and velocity field at 1200 Earth days after the simulation, where the first 200 Earth days 
were disregarded. Thus, the snapshots are of the exoplanet at about 34 HD 209458b days . The f our figure panels show the flow at different 
pressures and are chosen to approximately match Figures 1 and 2 of iRauscher & Menoul ( Eoifl) . who presented similar plots at P = 2.5 
mbar, 220 mbar, 4.4 bar and 20 bar at 1450 HD 209458b days. O ur plot for P = 2.13 mbar (top left panel) shows an upper atmosphere 
dominated by radiative forcing, similar to the top panel of Figure 1 of Rauscher & Menoul 1 201(3) . Farther down in the atmosphere at P = 216 



© 201 1 RAS, MNRAS 000. [71(231 



Atmospheric circulation of tidally -locked exoplanets 15 




1075 1139 1205 1267 1332 1396 1460 1525 1589 165317 1 8 1098 1155 1212 1269 1326 1383 1440 1497 1554 1611 1668 




100 200 300 100 200 300 

(deg) (deg) 



1062 1124 1186 1248 1310 1372 1434 1496 1558 1620 1683 1084 1142 1201 1260 1319 1378 1436 1495 1554 1613 1672 




100 200 300 100 200 300 

(deg) (deg) 



Figure 9. Snapshots of the flow field at about 340 HD 209458b days and at P = 216 mbar. Shown are results from the T21L33 (top left), T21L66 (top right), 
T31L33 (middle left), T31L66 (middle right), T63L33 (bottom left) and T63L66 (bottom right) spectral simulations. The hyperviscosity v is kept fixed such 
that the dissipation time varies and scales up when the resolution coarsens (see equation 1131 ). Colors indicate temperature in K, while arrows represent the 
velocity vectors. In these plots, the color bar range is fixed for clarity of comparison. 



mbar, a chevron-shaped fe ature is displaced eastwards of the substellar point and should be compared to the bottom panel of Figure 1 of 
Rauscher & Menou (2010). At P = 4.69 bar (bottom left panel), the advective time scales start to become shorter than r ra d, resulting in 



the longitudinal homogenization of temperature. At P = 21.9 bar (bottom right panel), the equatorial wind becomes more counter-rotating 
than super-rotating, partly as a result of the conservation of total angular momentum (which is set by starting the simulation from a windless 
initial state, in the absence of drag). Overall, th ere is a good degree of qua litative and quantitative agreement between our computed flow 
fields and those presented in Figures 1 and 2 of Rauscher & Menou j201fj) . despite the snapshots being taken at different times. However, 
some discrepancies remain: th e velocity features in our si mulations are stronger at the various vertical heights; the temperature ranges are 
discrepant from those shown in iRauscher & MenovJ (hold) , especially at P = 2.13 and 216 mbar. 
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Figure 10. Same as Figure|5] but for the finite difference core with K = 0.35. Shown are the G24L33 (top left), G24L66 (top right), G36L33 (middle left), 
G36L66 (middle right), G72L33 (bottom left) and G72L66 (bottom right) simulations. 



Figure [9] focuses on the P — 216 mbar snapshot at about 340 HD 209458b days, but for six different simulation resolutions: T21L33, 
T21L66, T31L33, T31L66, T63L33 and T63L66. We note that P ~ 0.1 bar is the pressure/height at which the infrared emission emerges and 
where the stratosphere (if any) begins. Details concerning the numerical resolution are given in Table[2] In general, the chevron- shaped flow 
feature is seen at all six resolutions — its substructure shows up in all of the simulations and is (expectedly) most clearly visible at T63. There 
are hints that the details of the flow, such as zonal wind speed, depend on the simulation resolution. The same conclusions may be drawn 
from the finite difference simulations presented in Figure [To] but we note that comparing B-grid simulations at different resolutions, with 
the same value of /C, may not constitute a fair exercise (see Appendix[C}. In other words, Figure|9]displays results from spectral simulations 
which are equally dissipative (same value of v) and at progressively finer resolutions, whereas Figure [Tolshows results from finite difference 
simulations which are more dissipative (numerically viscous) at lower resolutions. 

Other points deserve to be emphasized. Firstly, the spectral simulations manifestly capture the details of fine flow features better than 
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Figure 11. Minimum and maximum temperatures at P = 216 mbar, for both spectral (left panel) and finite difference (right panel) simulations, as functions 
of time in Earth days and as computed from the simulation snapshots taken in Figures l9land l 101 
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Figure 12. Temporally averaged, zonal-mean zonal wind (top row) and temperature (bottom row) profiles simulated for HD 209458b. Results from the spectral 
(left column; T63L33) and finite difference (right column; G72L33) simulations are shown. The horizontal dissipation parameters take their fiducial values of 
tv = 10 — 5 HD 209458b day and K, = 0.35. Temperatures are in K and wind speeds are in m s — 1 . 



the finite difference simulations, which is likely because the latter (with /C = 0.35) are more dissipative than the former (with t v = 10 
HD 209458b day). Secondly, there are clear qualitative differences between the snapshots from the spectral and finite difference simulations, 
which show ~ 10% variations in the temperature field. The third point concerns model variations and the level of variability in the temperature 
field at P = 216 mbar, which we illustrate in FigureQT] For the spectral simulations (left panel of Figure [Til, it is apparent that there are 
~ 10% variations in the temperature field as a function of time and the magnitudes of the variations are roughly equal across resolution (T21- 
T63). For the finite difference simulations (right panel of Figurell lb. the G72 simulations also show ~ 10% variations in the temperature field 
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Figure 13. Maximum values of the temporally averaged, zonal-mean zonal wind speeds, from different simulations, as functions of the vertical pressure, for 
various magnitudes of the horizontal dissipation. Here, "day" refers to one HD 209458b day which is about 3.5 Earth days. 
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Figure 14. Same as Figure [T2l but with the values of the horizontal dissipation parameters adjusted (t u = 10 6 HD 209458b day and K = 0.90) such that 
the zonal wind profiles (as shown in Figure [TSI approximately match. Temperatures are in K and wind speeds are in m s _1 . 



at P = 216 mbar whereas the G24 and G36 simulations show less variation, consistent with our earlier statement that the latter simulations 
are more dissipative (compared to their G72 counterpart). It is worth noting that numerical noise may contribute to these variations and thus 
complicate the comparisons. The level of variability present and its dependence on the magnitude of the horizontal dissipation applied will 
have implications for the study of variability in these model atmospheres for hot Jupiters. The discrepancies in the predicted temperature fields 
may be due to the different temporal evolutions within each simulation — even a pair of identical models seeded with slightly different (and 
random) initial perturbations may lead to different temporal evolutions. Furthermore, the time taken to reach quasi-equilibrium is different 
for each method of solution. 



© 201 1 RAS, MNRAS 000. [711231 



Atmospheric circulation of tidally -locked exoplanets 19 



Figure [T2l shows the usual Held-Suarez statistics for HD 209458b. The temporally-averaged, zonal-mean zonal wind profiles depict an 
equatorial, super-rotating wind down to ~ 10 bar, flanked by counter-rot ating jets at mid-latitude. W ind speeds have typical magnitudes 
~ 1 km s -1 . The top row of Figure [T2l shou ld be compared to F i gure 3 of iRauscher & Menoul (2010). While the temporally averaged and 
zonal-mean temperatures are not presented in lRauscher & Menoul 1 20101) . we still present these figures (bottom row) for comparison between 
the spectral and finite difference simulations. The temperature profiles between our pair of simulations are in good agreement. The similarity 
of the temperature profiles indicates that we have chosen ^discard (=200 Earth days) to be sufficiently large, such that differences due to 
initialization have been erased. However, noticeable differences exist particularly with respect to the wind field — the maximum speed of the 
equatorial, super-rotating wind is about 5 km s -1 in the finite difference simulation, but is only about 3.6 km s _1 in the spectral simulation. 

The discrepancies in the predictions for the zonal wind profile and the maximum speed of the super-rotating jet motivate us to explore 
the issue further by varying the values of t v and AC — the magnitude of horizontal dissipation (see Sj3 . 3 b — and subsequently varying 
the simulation resolution. Figure [T3] shows T63L33 and G72L33 simulations with various values of t v = 10~ r -10 -3 HD 209458b day 
and AC = 0.1-1, respectively. We note that AC = 1 is not used in normal circumstances, except near the poles to prevent the numerical 
problems previously described[f| The key point is that there are > 10% uncertainties associated with the predictions for the zonal wind 
profiles within each method of solution and also between them. The wind profile from the finite difference simulation with /C = 0.9 appears 
to approximately match that from the spectral simulation with t v = 10~ 6 HD 209458b day. Figure [T4l shows the Held-Suarez statistics 
with these adjusted values of the horizontal dissipation parameters — it is clear that there is now closer agreement between the pairs of 
simulations, both qualitatively and quantitatively. Overall, Figures Q3] and [14] demonstrate that one can, by trial and error, obtain consistent 
results for arbitrarily adjusted values oft v and AC, but there is still no rigorous way to choose these dissipation times/rates. 

To investigate the effects of varying the numerical resolution, we revert to the fiducial values of AC and the hyperviscosity v for the finite 
difference and spectral simulations, respectively. Keeping v constant while varying the resolution results in the variation of the numerical 
dissipation time/rate, thus allowing the spectral simulations to be compared on an equal footing (see equation 1131 ). However, keeping AC at a 
fixed value is strictly speaking analogous to varying v (see Appendix[C]l, so examining finite difference simulations at different resolutions, 
with the same value of AC, may not constitute a fair comparison. Nevertheless, we show in Figure [15] the results of our resolution studies 
from both the spectral and finite difference simulations, where we record the maximum speed of the temporally averaged, zonal-mean zonal 
wind as a function of P. For the spectral simulations, there is a spread of about 50% in the predicted wind speeds. The T31 and T63 
simulations agree reasonably well with one another at a given vertical resolution (L33 or L66), whereas the wind speed predictions from the 
T21 simulations are discrepant with the T31 and T63 ones, suggesting that the T21 simulations are under-resolved. For the finite difference 
simulations, the predictions for the wind speeds are also in agreement at a given vertical resolution (L33 versus L66). Simulations with lower 
horizontal resolution predict lower wind speeds, consistent with our earlier statement that they are more dissipative (at a fixed value of AC). 
Nevertheless, the depth at which the super-rotating wind is the fastest is a robust feature of the simulations, occurring at about 0.05 bar in the 
finite difference simulations (versus about 0.08 in the corresponding set of spectral simulations). 

We note that inter-comparison of the results between each panel of Figure [T5] does not constitute a fair exercise, because the assumed 
magnitude of horizontal dissipation is different for each suite of simulations. Even the intra-comparison of results within the right panel of 
Figure [15] may not be straightforward, because as we discussed previously keeping AC fixed while varying the resolution in effect changes 
the analogue of the hyperviscosity, but we do not know of a clear way of varying AC in a "resolution independent" manner (see Appendix 
O. The best we can conclude from Figure [T5] is that the predictions for the maximum zonal wind speed, from both the spectral and finite 
difference simulations, are resolution-dependent. In tandem with FigureQ~3] one may also conclude that since the specification of horizontal 
dissipation is a more lucid endeavour within the spectral core, comparing results from the spectral and finite difference simulations should 
only be performed when Held-Suarez statistics from the latter are calibrated to match those produced by the former. In this sense, the spectral 
simulations are more robust. 

We conclude that while we have achieved a rather satisfactory level of agreement between our spectral and finite difference simulations 
of HD 209458b, discrepancies arise in the quantitative predictions which may limit our ability to accurately model these extreme atmospheres, 
especially in terms of wind speeds. The main lessons we learn are that there are ~ 10% uncertainties associated with the temperature field 
and > 10% uncertainties associated with the velocity field, due to the choice of the magnitude of the horizontal dissipation as well as the 
resolution of the simulations. 



5 DISCUSSION 

5.1 Broader Implications 

Using a single and consistent simulation platform, we have performed a suite of benchmark tests concerning the atmospheric circulation 
of Earth and tidally-locked extrasolar planets. We find that while the dynamical cores of the FMS produce qualitative and quantitative 
agreement for the Earth, tidally-locked Earth and shallow hot Jupiter tests, the agreement is less than satisfactory for a deep model of the hot 
Jupiter HD 209458b. Further investigation reveals that closer agreement can be attained by arbitrarily adjusting the values of the horizontal 



We also examined a simulation with no horizontal mixing applied (AC = 0), which did not complete successfully (produced multiple output values of NaNs). 
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Figure 15. Maximum values of the temporally averaged, zonal-mean zonal wind speeds, from different simulations, as functions of the vertical pressure. Left: 
spectral simulations with different resolutions. Right: finite difference simulations with different resolutions. The fiducial value of the horizontal dissipation 
parameter used is K = 0.35 in the finite difference simulations, while the hyperviscosity v is kept fixed in the spectral simulations. Inter-comparing the left 
and right panels does not constitute a fair exercise, because the assumed magnitude of horizontal dissipation is different; even the intra-comparison of results 
within the right panel may not be straightforward (see text). 



dissipation parameters in the two dynamical cores, but there is no rigorous way to pick the magnitude of the horizontal dissipation in these 
models. 

Our findings suggest that even without dealing with additional physics such as radiative transfer or atmospheric chemistry, discrepancies 
in the temperature and velocity fields, at the level of 10% and several tens of percent respectively, al ready exist for the d ynamics alone. In 
this context, direct measurements of wind velocity in a hot Jupiter atmosphere, as recently reported by Snel len et alJ OOlOh . are important as 
they could prove to be particularly constraining for the models. 

In general, weakly-dissipative spectral simulations are expected to be sensitive to the choice of horizontal dissipation parameter (t v ), in 
the sense that they tend to fail as a result of small-scale noise accumulation ("spectral blocking") if the horizontal dissipation is not chosen to 
be strong enough. Therefore, to the extent that spectral simulations with the largest possible value of t v (i.e., the weakest possible horizontal 
dissipation) are more trustworthy — as conventional wisdom would suggest — our findings could also be interpreted as indicating that results 
in the literature based on finite difference methods may somewhat over-estimate the magnitude of wind speeds in h ot Jupiter atmosp heres 
(see Figure [T3l. However, until the nature of horizontal dissipation in these atmospheres is better understood (e.g.. lGoodmanl 2003), one 
should probably not interpret these trends as more than suggestive. 

Operationally, our suite of benchmark tests provides a reference for researchers wishing to adapt their codes to simulate atmospheric 
circulation on tidally-locked extrasolar planets, regardless of whether the codes solve the primitive or full Navier-Stokes equations. 



5.2 Summary 

The salient points of our study are: 

• We have generalized the Held-Suarez dynamical benchmark for Earth to include tidally-locked exoplanets using a single simulation 
platform (the FMS). Our suite of benchmark tests provides a reference for researchers wishing to adapt their codes to study atmospheric 
circulation on tidally-locked Earths/Neptunes/Jupiters. 

• We have found that the differences in the HD 209458b simulations of lCooper & Showman! j2005l.l200f3) and lRauscher & Menoul feoifi) 



are probably due to initial/boundary conditions and setup, and not due to the method of solution utilized. 

• Qualitative and quantitative agreement between the spectral and finite difference simulations of the deep- atmosphere benchmark test 
for the hot Jupiter HD 209458b can be attained if arbitrarily adjusted values of the horizontal dissipation parameters are adopted. However, 
the difficulty remains that the magnitude of the horizontal dissipation cannot (yet) be specified from first principles. This in turn leads to 
dynamical uncertainties at the level of > 10% which limit our ability to accurately model these atmospheres, especially with respect to wind 
velocities. Direct wind measurements from transit observations of extrasolar planets should thus be particularly constraining for the models. 
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APPENDIX A: IMPLEMENTING UNEVEN VERTICAL LEVELS 

The ability to perform dynamical core simulations with arbitrary spacing of the vertical levels comes as a default in the FMS, accom- 
plished using the following prescription: 

C = i-(*-i)/JV., 

C = A( + (1-A)C B , (Al) 
a — exp ( — QC J 



where the index i runs from 1 to N v + 1. Varying the parameters A, B and C allows one to control the range of a covered and the spacing 
between the points. For example, Fieure lATI shows three implementations of equation JAlb . where log a is evenly spaced only for B = 1. 



APPENDIX B: POLYNOMIAL FITS FOR THERMAL FORCING OF HD 209458B 



To aid the reader in reproducing our results, we provide convenient polynomial fits to r ra( j, Tnight an d ^day The radiative relaxation 
time is approximated by a 4th order polynomial fit, 



lot 



ITT ) 



5.4659686 + 1.4940124P + 0.66079196P^ + 0.16475329P J + 0.014241552P 4 

oo, 



P < 10 bar, 
otherwise, 



where P = log(P/l bar). The preceding fit is valid for 10 fibar ^ P ^ 8.5 bar. 
The night and day side temperatures are given by 



and 



^night 

1 K 



I'd ay 

1 K 



1388.2145 + 267.66586P - 215.53357P 2 + 61.814807P 3 + 135.68661P 4 + 2.0149044P J 



-40.907246P 
5529.7168 



19.015628P 7 - 3.8771634P 8 - 0.38413901P 9 - 0.015089084P 1 
6504P + 4142.7231P 2 - 936.23053P 3 + 87.120975P 4 , 



2149.9581 + 4.1395571P - 186.24851P 2 + 135.52524P' 4 + 106.20433P 4 - 35.851966P 5 
-50.022826P 6 - 18.462489P 7 - 3.3319965P 8 - 0.30295925P 9 - 0.011122316P 10 , 



5529.7168 - 6869.6504P + 4142.7231P 2 - 936.23053P 3 + 87.120975P 4 , 
respectively. Our fits presented in equations dB2t and jB3\ are valid for 1 /ibar ^ P ^ 3488 bar. 



P < 10 bar, 
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P < 10 bar, 
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Figure CI. E-folding time for initial grid-scale noise to decay away in the finite difference (B-grid) core of the FMS, tic, as a function of the horizontal 
mixing coefficient /C. The measurements are identical for the G24, G36 and G72 simulations. We show measurements (and fitting functions) from horizontal 
mixing schemes of both the second (squares) and fourth (triangles) order in wind, while retaining a fourth order temperature scheme, but present only results 
associated with the fourth order wind scheme in the paper. 



For completeness, we also provide a polynomial fit to the solid curve presented in Figure 1 of llro et al.l (2005): 
T 

= 1696.6986 + 132.23180P - 174.30459P 2 + 12.579612P 3 + 59.513639P 4 

1 K. 

+ 9.6706522P 5 - 4.1136048P 6 - 1.0632301P 7 + 0.064400203P 8 



(B4) 



+ 0.035974396P 9 + 0.0025740066P 10 . 



This fit is valid for 1 ^ibar < P < 3488 bar. 



APPENDIX C: E-FOLDING TIME FOR DECAY OF INITIAL GRID-SCALE NOISE (B-GRID CORE) 

A plausible way of determining the analogue of the numerical dissipation time t v (for the spectral core) in the case of the finite difference 
core — we shall denote this by tic — is to turn off all of the dynamical terms in the simulation except for damping and record the e-folding 
time for initial grid-scale noise (in velocity) to decay away. For completeness, we perform this task for a horizontal mixing scheme of both 
second and fourth order in wind, while retaining a fourth order temperature scheme at all times. Note that the simulation results presented in 
the paper are all associated with a fourth order scheme in both wind and temperature. 

Figure [CTl shows our measurements of tx. as a function of the horizontal mixing coefficient M As expected, the horizontal mixing 
scheme of second order in wind is more dissipative (smaller tic values) than the fourth order one. For convenience, we provide fitting 
functions to these results: 

t K _ J 3. 14235 - 5.22936/C + 9.61774£ 2 - 9.48662/C 3 + 3.55903/C 4 (2nd order wind), 

1 s [4.86355 - 10.2954/C + 19.3846£ 2 - 19.1606/C 3 + 7.21573/C 4 (4th order wind). ' 

The fits are valid for 0.1 ^ K, ^ 0.95 and only for the deep model of HD 209458b (cf. i]4.4l and Table[TJ. For example, with the wind scheme 
being fourth order, we get tic ~ 3 x 10~ 3 HD 209458b day for the fiducial value of K. = 0.35. 

It is important to note that the measurements of tic are identical for the G24, G36 and G72 simulations, implying that K, is strictly 
speaking the analogue of the dissipation time t v and not the hyperviscosity v, further implying that comparing spectral and finite difference 
simulations with different values of K, and t v (Figure [T3t is a meaningful exercise. Within the FMS, specifying K, constitutes a "resolution 
dependent" approach in that the smallest resolvable waves in the simulation are damped with a strength that is equal for every resolution. 
Unfortunately, there is no way within the FMS B-grid core (Memphis release) to switch to a "resolution independent" approach of specifying 
the magnitude of horizontal mixing, unlike in the case of the spectral core. It is also not clear how to convert K, into a quantity that is the true 
analogue of the hyperviscosity. The implication is that comparing simulations at different resolutions, with the same value of /C, may not be 
a fair exercise since the analogue of the hyperviscosity is different in each of these runs (Figure [10] and the right panel of Figurell5t. 

Even if we take Figure |CT1 at face value, the correspondence between t u and tic is not straightforward. For example, the pair of simu- 
lations in Figure [Uluse t v = 10~ 6 HD 209458b day and K. = 0.9, but FigurelCTl informs us that K. = 0.9 corresponds to t K ~ 4 x 10" 4 
HD 209458b day (using the horiziontal mixing scheme of fourth order in wind). If we insist on matching the K, = 0.9 simulation with the 
tu = 10~ 5 -10~ 4 HD 209458 day simulations, then Figure[T3linforms us that we are now faced with the original problem of the predictions 
for the maximum zonal wind speeds being discrepant. 
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